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Abstract This paper proposes a harmonic Lanczos bidiagonalization method for com- 
puting some interior singular triplets of large matrices. It is shown that the approxi- 
mate singular triplets are convergent if a certain Rayleigh quotient matrix is uniformly 
bounded and the approximate singular values are well separated. Combining with the 
implicit restarting technique, we develop an implicitly restarted harmonic Lanczos 
bidiagonalization algorithm and suggest a selection strategy of shifts. Numerical ex- 
periments show that one can use this algorithm to compute interior singular triplets 
efficiently. 
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1 Introduction 

The singular value decomposition (SVD) of a matrix A £ R MxN , M > N is given by 

A = UEV T , (1) 

where £ = diag{a\, • • • , (?n), U = (ui, u%, • • • , um) and V = (v±, V2, • ■ • , vjy) are 
orthogonal matrices of order M and N respectively. (<7j, Uj, «j), i = 1, 2, • • • N, are called 
the singular triplets of A. 
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Consider the (M + N) x (M + N) augmented matrix 

Then, the eigenvalues of A are ±01, ±02, • ■ ■ , ±crj\r and M — N zeros. The eigenvectors 
associated with cr^ and — cij are (itj , ^ J and (u^ , — v i J respectively. There- 
fore, the SVD problems are equivalent to the eigenproblems of augmented matrices. 

The SVD methods are widely used in determination of numerical rank, determi- 
nation of spectral condition number, least square problems, regression analysis, image 
processing, signal processing, pattern recognition, information retrieval, and so on. 

At present, computation of largest or smallest singular triplets of large matrices 
has been well studied, Lanczos bidiagonalization method and its variants are the most 
popular methods. In 1981, Golub et al. [5] firstly designed a block Lanczos bidiagonal- 
ization method to compute some largest singular triplets. Larsen [16] discussed the re- 
orthogonalization of the Lanczos bidiagonalization process. Jia and Niu [13] proposed a 
refined Lanczos bidiagonalization method to compute some largest and smallest singu- 
lar triplets. Kokiopoulou et al. [T5J used the harmonic projection technique to compute 
the smallest singular values. Baglama and Reichel [2][3] used Ritz values and harmonic 
Ritz values to approximate the largest and smallest singular values respectively. Her- 
nandez et al. [7] provided a parallel implementation of the Lanczos bidiagonalization 
method. Stoll [22] developed a Krylov-schur approch to partial SVD. Recently, Jia and 
Niu [14] proposed a refined harmonic Lanczos bidiagonalization method to compute 
some smallest singular triplets. All of above methods compute the Lanczos bidiagonal- 
ization process, build two m— dimensional Krylov subspaces, then extract approximate 
singular triplets from these two subspace by different ways. Hochstenbach 8,9 also 
give the Jacobi-Davidson type algorithms for SVD problems. 

Due to the storage requirement and the computational cost, all the projection meth- 
ods must be restarted. The implicit restarting technique 21 proposed by Sorensen is 
the most powerful tool and is widely used in many projection methods. The success of 
this technique heavily depends on the selection of the shifts, see [101121] . For eigenvalue 
problems, Sorensen [21] used the unwanted Ritz values as the shifts to restart Arnoldi 
method, and Morgan [19] used the unwanted harmonic Ritz values as the shifts to 
restart harmonic Arnoldi method. Jia [101111] used the refined shifts and refined har- 
monic shifts obtained by the information of the refined Ritz vectors and refined har- 
monic vectors to restart refined Arnoldi method and refined harmonic Arnoldi method, 
respectively. For SVD problems, Kokiopoulou et al. [T5J used the unwanted harmonic 
Ritz values as the shifts. Baglama and Reichel [2][3] explicitly augmented the Lanczos 
bidiagonalization method with certain Ritz vectors or harmonic Ritz vectors. Jia and 
Niu [131114] gave an refined (harmonic) shift strategy within the implicitly restarted 
refined (harmonic) Lanczos bidiagonalization method. 

In this paper, we are concerned with the computation of interior singular triplets. 
For a given target r, we want to compute some singular triplets nearest r. So, we sort 
the singular triplets by 

|<T1 — t| < |<72 — t\ < ■ • • < |<7jV — t\. (3) 

We must emphasize that, in this paper, a\ is the singular value nearest r rather 
than the smallest singular value, meanwhile, ctjv is the singular value farthest from r 
rather than the largest singular value. 
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Since the largest eigenvalues of (A — tI)^ 1 are the eigenvalues of A closest to r, and 
the SVD problem of A is equivalent to the eigenproblem of A, we can use shift-invert 
technique on A — rl to compute the interior singular triplets, such as shift-and-invert 
Arnoldi method (svds). In this paper, we assume that M and N are large and that A 
can not be factorized. The shift-and-invert technique need the factorization of A — rl. 
Since M and N are large, M + N, the dimension of A — tI, is larger. We can not do 
any factorizations on A — rl. Therefore, the shift-and-invert technique is not suitable 
for interior SVD problems. 

Another approach for computing interior singular triplets is the harmonic pro- 
jection method. The harmonic projection method has been widely used to compute 
interior eigenpairs, see [181119] . and has been adopted to combine with Lanczos bidiag- 
onalization methods to compute smallest singular triplets [2ll3lll5lll3| . However, if we 
use the harmonic projection method explicitly on A — tI, the scale of the problem is 
increased and this leads to the increasing computational cost. Further, we ignore the 
special structure of A or A — rl, and the projected matrix and the updated process 
of implicit restarting may lose this structure. Therefore, we must use the harmonic 
projection method implicitly. Until now, no literature has been appeared to compute 
interior singular triplets by the harmonic projection method implicitly. 

In this paper, we propose a harmonic Lanczos bidiagonalization method for com- 
puting interior singular triplets by combining the harmonic projection technique with 
the Lanczos bidiagonalization process. We analyze the convergence behavior, show that 
the harmonic Ritz approximations converge to the desired interior singular triplets if 
some Rayleigh quotient matrix is uniformly bounded and the harmonic Ritz values are 
well separated. Then, based on Morgan's harmonic shift strategy [TU] for computing 
interior eigenvalues, we give a selection of the shifts within the framework of the implic- 
itly restarted harmonic Lanczos bidiagonalization methods. Further, we report some 
numerical experiments of computation of interior singular triplets. It appears that the 
algorithm we proposed is suitable for computing the interior singular triplets of large 
matrices. 

Throughout this paper, denote by || ■ | the spectral norm of a matrix and the vec- 
tor 2-norm, by K. m {C,vi) = span{vi, Cv\, ■ ■ ■ , C m ~ 1 vi} the m— dimensional Krylov 
subspace generated by the matrix C and the starting vector v\, by superscript ' T ' the 
transpose of matrix or vector, by e m the m— th coordinate vector of dimension m. 



2 Harmonic Lanczos bidiagonalization method 

2.1 Lanczos bidiagonalization process 

Golub et al. [5] proposed a Lanczos bidiagonalization method to compute the largest 
singular triplets of A. This method is equivalent to the symmetric Lanczos method on 
A with a special initial vector. It is based on the Lanczos bidiagonalization process, 
which is shown in matrix form as follows: 

AQm = P m B m , (4) 

^4 Pm — QmB m + /3 m 9m+l e mi (5) 



1 



where 



/ai Pi 



B m = 



«2 



V 



(6) 



is an upper bidiagonal matrix 



21,92, 



am ) 

,q m ) and P m = (pi,P2,- 



p m J span 

the Krylov subspaces /C m ( J 4 1 A, qi) and K, m {AA l ,pi), respectively. 

In finite precision arithmetic, the columns of P m and Q m may lose the orthogonality 
rapidly and must be reorthognalized. From the analysis of Simon and Zha [2D], we know 
that only the columns of one of the matrices P m and Q m need to be reorthogonalized. 
When M 3> N, Reorthogonalization on Q m only can reduce the computational cost 
considerably. So we only perform reorthogonalization on Q m - 



2.2 Harmonic Lanczos bidiagonalization method 
Given the subspace 

£ = ^"{(tol)}- < r > 

Making use of the harmonic projection principle, we compute some approximate eigen- 
pairs {6i,<f>i) of A nearest r by requiring 

(A-eiI)<pi±(A-Tl)S. 1 ' 

From © and (O, © can be rewritten as the following generalized eigenproblem: 

2 t i d dT i o „ „T 



r/ B m \ /iA 1 f t I + B m B m + Pme m e m -2rB m \ ( x 



Bm -rl ) \Vi J 6i-r\ -2rB m T 2 I + B m B m J \yi 



(9) 



Assume that Bi > 0, i = 1, 2, • • • , k + Z, which are sorted by 
|6»i-r|<|6» 2 -r|<---<|6> fc+i -T| 

and #i < 0, i = fe+i+1, k+l+2, • ■ ■ , 2m. We can use i = 1, 2, • • • , k and = ( g m ^ ^ 

as the approximation of the desired eigenpair of A. Because of the relation between the 
singular triplets of A and the eigenpairs of A, we use Qi,U{ = P m Xij\\xi\ \ — P m Xi,Vi = 
Qmyi/\\yi\\ = QmVi, i = 1,2, •• •, k as the approximate singular triplets of A nearest r. 
Here we call B{,ui, Vi the harmonic Ritz value, the left and right harmonic Ritz vector, 
respectively. 

From Q and (|5j), we have 

\\Avi - 6iUi\\ = \\BmVi - 0iXi\\, 



\\Avi- 0iUi\\ = \Zl\BmXi - OiViW 2 + Pm\e m Xi\ 



I - ■" j " : " i 1 1 — v 1 1 ■■' ' I 1 " ' '" 

Therefore, if 



y/\\BmVi - ^|| 2 + \\B&Xi - GmW 2 + I3 m \e7n£ l \ 2 < tol, (10) 

where tol is a prescribed tolerance, then the method is known as convergent. So we 
need not form -Sj and 5j explicitly before convergence. 
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2.3 Convergence analysis 
Set 

B = 

and 



-rl B m 
Bl -rl 



C 



T 2 I + B m Bm + /3 m e m em -2TB m 

-2tB^ T 2 I + Bj n B m 



then 6i,i = 1,2, •••,2m are the eigenvalues of B~~ l C. The matrix B is called the 
Rayleigh quotient matrix of A with respect to the subspace £ and the target r. 

The following results are direct from Theorem 2.1, Corollary 2.2 and Theorem 3.2 
of [12]. 

Theorem 1 Assume that (a, u, v) is a singular triplet of A, define that e = sin Z ((") ,£ ) 
is the distance between the vector ("j and the subspace £. Then there exists a pertur- 
bation matrix F such that a is an exact eigenvalue of B~ l C + F, where 

llFll^^^jjB^IKallAII + IIAH 2 ). (11) 

Furthermore, there exists an eigenvalue of B~^C satisfying 

|0-<7|<(2||A|| + ||F||)||F||. (12) 

Theorem [T] shows that if e tends to zero and if ||_B _1 || is uniformly bounded, then 
there exists one harmonic Ritz value 9 converging to the desired singular value a. 
However, from the interlacing theorem of eigenvalues [6], since 

5 ( Pm \ T - (P m 

B= { Q m ) ( A - Tl \ Q m 

we have that the eigenvalues of B are between the largest and smallest eigenvalue of 
A — tI. Therefore, B may be singular, which leads to arbitrarily large ||_B _1 ||. Hence, 
we must assume ||_B _1 || is uniformly bounded. In fact, this is the inherent defect of the 
harmonic projection methods, which can be easily obtained from Jia's analysis [12j . 

Similarly to the analysis in [T2], if r is very close to a desired singular value a of A, 
then the method may miss it. We replace 6^ by the Rayleigh-quotient pi = uj Avi = 
xj Brajji as the approximate singular value, as was done in [9lll4], In general, pi is more 
accurate than 9i. 

Theorem 2 Let (9,z) be an eigenpair of B~ X C, where z = (yj> an ^ assume i z ,Zj_) 
to be orthogonal such that 



zl) B c ^ = [og)- (13) 

If 

sep(8,G) = \\{G-9ir 1 \\- 1 >0, (14) 
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then 



Vl — e 2 sep(a, G) 

<(i + - m V- (15) 

Theorem [5] shows that if ||_B _1 || is uniformly bounded and sep(9,G) is bounded 
below by a positive constant, that is, all harmonic Ritz values are well separated, then 
the harmonic Ritz vectors u, v converge to the desired left and right singular vector. 



3 Implicit restarting technique, shifts selection and an adaptive shifting 
strategy 

3.1 Implicit restarting technique 

Due to the storage requirement and the computational cost, the number of Lanczos 
bidiagonalization steps m can not be large. However, for a relatively small m, the 
approximate singular triplets do not converge. Therefore, the method must be restarted 
generally. 

The implicit restarting technique proposed by Sorensen |21j is a powerful restarting 
tool for the Lanczos and Arnoldi process, and has been adopted to the Lanczos bidi- 
agonalization process [4 , 13, 14, 15, 17 . After running the implicit QR iteration p steps 
on B m and using the shifts fij, j — 1, 2, • • • ,p, we have 

{B m B m — jttfi) • ■ • (B m B m — (J,pl) — PR, 

P^B-mQ upper bidiagonal, 

where P, Q are the products of the left and right Givens rotation matrices applied to 
B m - 

Performing the above process gives the following relation: 

AQm—p Pm—p^m—pi (17) 
A P m - P = Qm-pBm—p + {Pm—pPm,m-p1m+l + Pm—p'lm-p+l) e ™-p-> 

where Q m _ p and q m _ p+ i are the first m — p columns and the (m — p+ l)-th column of 
QmQ, Pm-p is the first m — p columns of P m P, B m _ p is the leading (m — p) x (m — p) 
block of PB m Q, p m ,m- P is the (m,m — p) element of P. Since P m -pPm,m-p<lm+X + 
Pra— p%n— p+1 ls or th°g ona l to Q m _ p , we obtain a (m — p)-step Lanczos bidiagonaliza- 
tion process starting with , where 

p 

lit = Y[(A T A-^I) qi (19) 
i=i 

with 7 a factor making ||g^|| = 1. It is then extended to the m-step Lanczos bidiago- 
nalization process in a standard way. 
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3.2 shifts selection and adaptive shifting strategy 

Once the shifts (4%, ji2, ■ ■ ■ , [ip are given, we can run the implicitly restarted algorithm 
described above iteratively. The success of the implicit restarting technique heavily 
depends on the selection of the shifts. As is shown in [13], from ()19[) . it can be easily seen 
that the more accurate the shifts approximate to some unwanted singular values, the 
more information on the unwanted singular vectors are dampened out after restarting. 
Therefore, the resulting subspace contains more information on the desired singular 
vectors, and the algorithms may converge faster. For eigenproblems and SVD problems, 
Morgan [19] and Kokiopoulou et al. [15] suggested using the unwanted harmonic Ritz 
values as shifts. A natural choice of the shifts within our algorithm is the unwanted 
approximate singular values 6^ + j ,j = 1, 2, • • • , I, since they are the best approximations 
available to some of the unwanted singular values within our framework. 

From (|19p , we see the component along the desired fe-th singular vector uj, is greatly 
damped if a shift fij is very close to a^, so fij is a bad shift and p^ may converge to 
<jj- very slowly or not at all. To correct this problem, Larsen |17j proposed an adaptive 
strategy to compute largest singular triplets. He replaces a bad shift by zero shift. Jia 
and Niu 13 II gave a modified form for computing smallest singular triplets. Define 
the relative gaps of p^ and all the shifts pi, i = 1, 2, • • • , I by 



relgapj 



{Pk - £ k) - Mi 



Pk 



(20) 



where is the residual norm (|10|l . If relgap fci < 10~ 3 , m is a bad shift and should 
be replaced by a suitable quantity. They replace the bad shifts by the largest or the 
smallest approximate singular value for computing the smallest or the largest singular 
triplets. In this paper, a good strategy is replacing the bad shifts by the approximate 
singular value farthest from r, as this strategy amplifies the components of in 
Vi,i = 1, 2, • • • , k and damps those in Vi,i = k + 1, k + 2, • • • , N. 



4 Numerical Experiments 

Numerical experiments are carried out using Matlab 7.1 R14 on an Intel Core 2 E6320 
with CPU 1.86GHZ and 2GB of memory under the Window XP operating system. 
Machine epsilon is e mac h ~ 2.22 x 10~ 16 . The stopping criteria is 



If 



stopcrit = max yjjACj — PiUi\\ 2 + || J 4 T fii — piVi\\ 2 . (21) 

l<i<k 

w < tol > (22) 

then stop. From (|10[1 . we need not form Ui,Vi explicitly before convergence. 

For large eigenproblems, in order to speed up convergence, most of the implicitly 
restarted Krylov type subspace algorithms, such as ARPACK(eigs), compute k + 3 
approximate eigenpairs when k eigenpairs are desired. This strategy has been adopted 
to SVD problems, see [2lll4|. In this paper, we also compute A: + 3 approximate singular 
triplets and use / — 3 shifts in implicit restarting process. 

All test matrices are from pQ. We take tol = 10 -6 . In all the tables, ''iter'' denotes 
the number of restart, 'time' denotes the CPU timings in second, 'mv' denotes the 
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Fig. 1 Absolute residual norms for WELL1850 for k = 3, m = 20, r = 0, 0.01, 0.005, 0.001 



number of matrix-vector products. Since the matrix-vector products performed on A 
are equal to those on A T , we only count the matrix- vector products on A. 



4.1 Computation of smallest singular triplets 



Obviously, we can compute some smallest singular triplets by taking r = 0. We com- 
pute three singular triplets nearest r = 0, 0.01, 0.005, 0.001 of WELL1850, respectively. 
These three singular values are all the three smallest singular values. The computed 
three singular values are 



a x « 1.611969e - 002, a 2 « 1.911309e - 002, ct 3 « 2.315889e - 002. 



Table [T] reports the computational results. Fig. [T] plots the absolute residual norms 
of the computed singular triplets for m = 15 and m = 20, respectively. From Table 
Q] and Fig. [T] we see that for all r, our algorithm can compute three singular triplets 
accurately. However, for different r, the algorithm has a great difference on restart 
numbers, matrix-vector products and CPU times. This phenomenon shows a good 
choice of target point r can speed up the convergence considerably. 



Table 1 WELL1850 for k = 3, m = 10, 15, 20, 25, r = 0, 0.001, 0.005, 0.01 



m 


iter | time \ mv \ stopcrit 


iter | time \ mv \ stopcrit 




T — 


t = 0.001 


10 


543 


5.15 


2178 


1.67c-005 


178 


1.89 


718 


1.66e-005 


15 


119 


3.20 


1077 


1.67c-005 


74 


1.99 


672 


1.25e-005 


20 


56 


3.48 


790 


1.50o-005 


48 


2.81 


678 


1.62e-005 


25 


35 


3.20 


671 


1.35o-005 


35 


3.64 


671 


1.14e-005 




t = 0.005 


T = 0.01 


10 


160 


1.66 


646 


1.66c-005 


179 


1.76 


722 


1.60c-005 


15 


68 


1.86 


618 


1.62c-005 


64 


1.73 


582 


1.52c-005 


20 


39 


2.31 


552 


1.08c-005 


37 


2.23 


524 


1.45c-005 


25 


31 


3.10 


595 


1.27c-005 


29 


2.89 


557 


1.06c-005 
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4.2 Computation of three interior singular triplets nearest different r 

The test matrix is DW2048, a 2048 x 2048 matrix. We compute three singular triplets 
nearest different r. The computational results are shown in Tables [2][3] From Table 
[21 we see that the relative errors of the computed singular values are no more than 
O(10 -9 ). The Tables demonstrate that our algorithm can compute the desired singular 
triplets accurately. 



Table 2 Three computed singular values of DW2048 nearest r = 0.2, 0.5, 0.6, 0.8 for m = 50 



t = 0.2 


r = 0.5 


Pi 


\pj - CTjI/CTj 


Pi 


\pj - CTj|/CTj 


2.0031301c-001 


1.55e-014 


4.9933773c-001 


1.62c-14 


1.9939880c-001 


5.90c-014 


5.0082218c-001 


1.04e-12 


1.9813769c-001 


1.08c-009 


4.9764898c-001 


8. 76c- 11 


t = 0.6 


t = 0.8 


Pi 


\pj - (Tj\/(Tj 


Pi 


\PJ - CT l\/ CJ 3 


6.0106012c-001 


4.29e-12 


8.0014466c-001 


2.41c-12 


6.0193472e-001 


4.22e-ll 


7.9954438c-001 


5.46e-12 


5.9689466c-001 


2.40e-13 


7.9932106c-001 


1.08c-10 



Table 3 DW2048 for k = 3, m = 30, 40, 50, r = 0.2, 0.5, 0.6, 0.8 





t = 0.2 


t = 0.5 


m 


iter 


time 


mv 


stopcrit 


iter 


time 


mv 


stopcrit 


30 


501 


109 


11655 


9.97c-007 


255 


51.2 


6123 


9.93c-007 


40 


298 


113 


9993 


9.81c-007 


97 


38.2 


3271 


9.90e-007 


50 


221 


136 


9652 


9.85c-007 


83 


50.9 


3656 


9.10c-007 




r = 0.6 


r = 0.8 


m 


iter 


time 


mv 


stopcrit 


iter 


time 


mv 


stopcrit 


30 


125 


24.6 


3006 


9.11e-007 


405 


78.2 


9525 


9.87c-007 


40 


69 


27.0 


2343 


9.47c-007 


180 


70.4 


6079 


9.64c-007 


50 


46 


28.1 


2012 


8.61e-007 


136 


81.8 


5989 


9.50e-007 



4.3 Computation of interior singular triplets for different k 

We compute k — 1,3,5, 10 smallest singular triplets nearest r = 4.5 of LSHP2233, a 
2233 x 2233 matrix. Table [5]reports the results. We see that our algorithm can compute 
the desired singular triplets with high precision. 



5 Conclusion 

In this paper, combining the harmonic projection principle with the implicit restarting 
technique, we propose an implicitly restarted harmonic Lanczos bidiagonalization algo- 
rithm for computing some interior singular triplets. Based on Morgan's harmonic shift 
strategy for computing interior eigenpairs, we give a selection of the shifts within our 



10 



Table 4 Ten computed singular values of LSHP2233 nearest t = 4.5 for m = 50 



Pi 


Pi — <Tl|/<Tl 


P2 


P2 — 02I/02 


4.4988631 


1.58c-15 


4.5091282 


1.36e-14 


P3 


P3 —0-3/(73 


P4 


P4 - 04I/04 


4.5113859 


1.22e-14 


4.4815289 


6.54e-15 


P5 


P5 - 05 | / (T5 


P6 


P6 - o 6 |/o 6 


4.5188882 


5.11c-15 


4.5210494 


1.18e-14 


P7 


P7 — 0-7/0-7 


P8 


P8 - o s |/o 8 


4.4783693 


1.07e-14 


4.4716358 


8.74e-15 


P9 


P9 —09/09 


P10 


|pio — 010 /010 


4.5331457 


5.68e-15 


4.4638926 


2.03e-10 



Table 5 LSHP2233 for k = 1, 3, 5, 10, m = 30, 40, 50, r = 4.5 





fc= 1 


k = 3 


m 


iter 


time 


mv 


stopcrit 


iter 


time 


mv 


stopcrit 


30 


467 


108 


11920 


6.88e-006 


560 


127 


13171 


6.95e-006 


40 


190 


83.6 


6844 


6.86e-006 


230 


98.1 


7826 


6.79c-006 


50 


159 


114 


7188 


6.63e-006 


216 


140 


9404 


6.66e-006 




= 5 


k = 10 


m 


iter 


time 


mv 


stopcrit 


iter 


time 


mv 


stopcrit 


30 


322 


64.8 


6972 


6.99c-006 


651 


103 


10761 


6.91c-006 


40 


207 


78.5 


6632 


6.78e-006 


168 


61.0 


4548 


4.70c-006 


50 


132 


83.8 


5487 


6.40e-006 


165 


92.4 


6003 


6.99c-006 



algorithm. Numerical experiments show that our algorithm is suitable for interior SVD 
problems. The interior singular values can be computed with higher relative precision. 
The Matlab code can be obtained from the authors upon request. 

Acknowledgements We thank Baglama and Reichel very much for generously providing 
their Matlab code of Lanczos Bidiagonalization process on their homepage, which reduces our 
programming work greatly. 
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